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ID and 2D simulations for the NASA Electric Arc 
Shock Tube experiments 

By D.V. Kotov, H.C. Yee, M. Panesi, D. Prabhu, and A. Wray 


1. Motivation and objectives 

The Electric Arc Shock Tube (EAST) facility at NASA Ames Research Center is used 
to generate high-enthalpy gas tests for studying high-speed atmospheric entry physics. 
The facility is composed of a long tube and a chamber. In the chamber, which contains 
a driver gas, an electrical discharge takes place, giving rise to a sudden increase of the 
gas temperature and pressure. The chamber’s geometry can be approximated as a 10- 
cm cylindrical tube though its real geometry is more complex and not axisymmetric. 
An aluminum diaphragm separates the driver from the driven gas. At the high pressure 
generated by the discharge, the diaphragm bursts, forming a shock wave that travels at 
high speed through a long cylindrical tube, 10 cm in diameter. As the shock propagates 
downstream, the shock-heated gas radiates, and in a test section emission spectroscopy is 
used to determine the radiative signature and thereby the thermo-chemical and radiative 
properties of the medium. 

The experiments being simulated here make use of helium as a driver gas and synthetic 
air (N 2 + O 2 ) as a test gas (or driven gas). The shock velocities obtained in EAST 
experiments range between 9 and 16 km/sec. The distance between the diaphragm and 
the test section is 7.0 m. At the test section the spectrally and spatially resolved shock- 
layer radiance is analyzed by taking a snapshot of the shock wave and the following gas 
as they pass in front of an optical access window. 

It is important to note that optically probing the shocked gas does not provide any 
information about the radial structure of the flow in the shock tube. The experiments only 
provide integrated measurements that include all the absorption and emission across the 
tube. Experimental data processing currently conducted at NASA Ames takes the flow 
to be one-dimensional, and any boundary layer effects are neglected. In reality, however, 
viscous effects may significantly impact the interpretation of the tests. To make better 
use of the flow radiation measurements one needs to know the radiative properties within 
the boundary layer, which can be estimated by numerical simulations. 

However there are some challenges in a CFD computation of this problem. One major 
challenge lies in the specification of the thermochemical state of the arc-heated driver gas, 
i.e., the initial state. Although energy deposition into the driver gas is time dependent, 
the current approach has been to adjust a constant temperature until the target shock 
velocity is reproduced. 

Several other challenges arise from the need to accurately capture the spatial-temporal 
evolution of the flow field in order to properly characterize boundary layer growth on the 
shock tube wall. The boundary layer is the primary cause of the shock wave’s decelera- 
tion, and it has important implications for re-absorption of shock layer radiation prior 
to instrument detection. These phenomena are unsteady and must be treated as such. 
For flight experiments one usually seeks a steady state solution, so that the simulation 
is not constrained by time accuracy requirements, and various acceleration techniques 
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can be used. The shock tube flowfield, however, develops in time and thus requires, at 
a minimum, second order time integration and minimization of numerical instabilities 
due to the multiscale physics and stiff source terms of the flow. Similarly, the spatial 
relaxation scale has to be resolved accurately, and all wave interactions have to be ac- 
counted for. Initial simulations of the EAST facility have demonstrated strong sensitivity 
to the choice of numerical method. Therefore, schemes with high-order spatial accuracy 
and low dissipation are desired. Furthermore, the length of the EAST facility is roughly 
8 m, which translates to enormous grids. The combination of these temporal and spatial 
resolution requirements demands significant computer resources and simulation time to 
account for all important phenomena. 

An additional, less well-known challenge involves the stiffness of the source terms in 
the governing equations. The temperature in the shock during the experiment is above 
20, 000 K, which means that chemical reaction rates in the vicinity of the shock be- 
come very high and the source term describing the chemical reactions becomes stiff. As 
most common shock-capturing schemes have been developed for problems without source 
terms, when applied to problems with nonlinear and/or stiff source terms these methods 
can result in spurious solutions, even when solving a conservative system of equations 
with a conservative scheme. This kind of behavior can be observed even for a scalar case 
(LeVeque & Yee 1990) as well as for the case with two species and one reaction (Wang 
et al. 2012). For further information concerning this issue see (LeVeque & Yee 1990; 
Griffiths et al. 1992; Lafon & Yee 1996; Yee et al. 2011, 2012). 

For a brief introduction and earlier CFD simulations of EAST see (Yee et al. 2012; 
McCorkle & Hassan 2010). Because of the non-symmetrical geometry of the driver zone, 
partially due to the opening of the diaphragm and the non-symmetrical configuration 
of the discharge, a full 3D computation should be carried out for this problem. The 
present investigation is to perform less CPU-time-intensive ID and 2D computations for 
the purpose of gaining first-hand understanding of the simulation challenges involved. 
Knowledge gained will provide guidance for future 3D CFD simulations. 


2. Governing equations 

Consider the 3D reactive Navier-Stokes equations for a one-temperature energy model: 

dp d 

~dt ^ 'dx~-^ sU ^ ~ ^ s (2.1) 

d d 

^ (frui) + —(puiUj + pSij - nj) = 0 ( 2 . 2 ) 

d d 

-Q t (pE) + — (uj(pE + p) + Qj + y ^ p s d sj h s - UiTij) = 0, (2.3) 
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where i = 1,2,3, (p s , put, pE) are the conserved variables, p s are the partial densities 
with k = 1, . . . , N s for a mixture of N s species. The mixture total density, the pressure, 
and the total energy per unit volume are 

N s N s 

P = J2 Ps ' p= P E =^ J Ps{ e s(T) + h° s ) + -pv 2 , (2.4) 

s s=l 1 s S=1 

where R is the universal gas constant, h° s are the species formation enthalpies, and M s 
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are the species molar masses. The viscous stress tensor is 


Tij = p 


dui diij 

— l - 

dxj dxi 


2 du k 

M 3dx k Sij ' 


The diffusion flux is 


j n dX s 
83 ~ s d Xi 


(2.5) 


( 2 . 6 ) 


where D s is the diffusion coefficient and X s is the mole fraction of species s. The con- 
ductive heat flux is 

dT 

« = ^ {2n) 
where A is the thermal conductivity of the mixture. The chemical source term is 


N r 


£l s — -M s ^ ^ (ps,r &s,r) 


N s 
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N s 


Pn 

M„ 


( 2 . 8 ) 


where a and b are the stoichiometric coefficients, and the forward reaction rate coefficients 
kf r are given by Arrhenius’ law: 

k fr = A f ^T ntr exp (-Ef^r/kT) (2.9) 

The backward reaction rate coefficients are computed as kb, r = kf^/Kff. where Kf q r is 
the equilibrium constant 


3. Computation Results 

All the computations employ the multi-dimensional high order single/overset grid 
nonequilibrium code ADPDIS3D (Lani et al. 2013; Sjogreen & Yee 2009). ID and 2D 
simplifications of the 3D EAST problem are considered using a single block option in the 
code. In addition, a one-temperature model is used, though the two-temperature model 
is required in order to obtain better agreement with the experiments. For the 2D case, 
axisymmetric geometry is not used as this option has not been implemented into AD- 
PDIS3D; planar geometry is used instead. The MUTATION library (Panesi et al. 2011), 
developed by Thierry Magin and Marco Panesi, is used to provide reaction rates and 
transport properties. We shall compare the results obtained by several standard shock- 
capturing methods and their filter counterpart schemes (Yee & Sjogreen 2007, 2010) for 
the early time evolution of the flow. Note that, for this viscous simulation, all the CFL 
values are based on the convection and viscous parts of the PDEs. Unless indicated, all 
shock-capturing schemes use the Roe average states. 

3.1. ID EAST simulation results 

The computational domain has a total length of 8.5 to. The left side of the domain, with 
length 0.1 in, is a high pressure region. The right side of the domain, with length 8.4 to, 
is a low pressure region. The temperature in the vicinity of the shock can reach more 
than 20, 000 K , therefore ionized species must be taken into account. Here we consider 
the gas mixture as consisting of 13 species: e~ , He, N, O, N 2 , NO, 0 2 , N 2 , NO + , N + , 
0+ 0+, He+. 

The initial conditions of the high and low pressure regions are listed in the table 1. The 
initial driver gas temperature is taken to be 6000 K, as in McCorkle & Hassan (2010), 
and the pressure p = 12.7116 MPa is chosen to obtain a shock velocity of approximately 
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Figure 1. ID, 13 species EAST problem: Second-order Harten-Yee TVD simulation for four 
levels of grid refinement with CFL = 0.8 and t en d = 3.25 x 10 -5 sec. The Aa; grid refinement 
study is: 10 -3 m (line 1), 5 x 10 -4 m (line 2), 5 x 10 -5 m (line 3), 2.5 x 10 -5 m (line 4). 



Figure 2. ID, 13 species EAST problem: Comparison among 5 methods using 501 point grids 
with CFL = 0.8 and t en d = 3.25 x 10 -5 sec. Reference solution (TVD on a 10,001 point grid) 
(line 1), TVD (line 2), TVDafi+split (line 3), WEN05-llf (line 4), WEN05Pafi+split (line 5), 
TVD/SR (line 6). See text for method notation. 
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p 

1.10546 kg/m 3 

P 

3.0964 x 10 ~ i kg/m' 

T 

6000 K 

T 

300 K 

P 

12.7116 MPa 

P 

26.771 Pa 

Y He 

0.9856 

Yo 2 

0.21 

Yn 2 

0.0144 

Yn 2 

0.79 


Table 1. High (left) and low (right) pressure region initial data 



Figure 3. ID, 13 species EAST problem: shock velocity time dependance solution using 
At = 5 x 10 “ 'sec on four Ax: 10 -3 m (line 1), 5 x 10 -4 m (line 2), 5 x 10 -5 m (line 3), 
2.5 x 10~ 5 m (line 4). 


10 km/sec in our computation, which is a typical velocity observed in the experiment. 
For the left-side boundary an Euler (slip) wall condition is applied, and for the right-side 
a zero gradient condition is applied for all variables. 

Uniform 1-D grids are used for these simulations. To save computational cost, an initial 
computational domain of (—0.1, 0.4) m is generated. During the computation the shock 
location is calculated at each time step. When the shock is close enough to the the right 
boundary, the computational domain is increased on the downstream side by 0.5 m. 

Figure 1 shows the results from the computation using the Harten-Yee second-order 
TVD scheme (Yee 1989; Yee et al. 1990) for four grids with Ax = 10 -3 m, 5x 10 -4 to, 5x 
10 -5 m and 2.5 x 10 -5 m at time t en( i = 0.325 x 10 -4 sec. One can observe a significant 
shift in the shear (left discontinuity) and the shock (right discontinuity) locations as the 
grid is refined. The distance between the shear and the shock shrinks as the grid is refined. 
The difference between shock locations obtained on the grids with Ax = 5 x 10 -5 to and 
2.5 x 10 -5 to is less than 0.3%. Thus the solution using Ax = 5 x 10 -5 to can be considered 
as the reference solution. 

Figure 2 shows a comparison between each of the five methods obtained on a coarse 
grid (Ax = 10 -3 ) to, along with the reference solution. The scheme labels are defined as 
follows: 
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Figure 4. Computational domain for the 2D EAST problem. 

• TVDafi+split: Sixth-order central base scheme with the Ducros et al. (2000) splitting 
of the governing equations. The flow sensor for the filter step is based on the shock and 
shear locations instead of using wavelets. See Yee & Sjogreen (2007, 2010) for further 
information on filter schemes. 

• WEN05-llf: Fifth-order WENO (WEN05) using the local Lax-Friedrichs flux. 

• WEN05Pafi+split: Nonlinear filter counterpart of the positive WEN05 using the 
local Lax-Friedrichs flux. The flow sensor for the filter step is based on the shock and 
shear locations instead of using wavelets. For the finite difference form of the positive 
WEN05-llf, see Zhang & Shu (2012). 

• TVD/SR: Finite difference scheme with subcell resolution (Wang et al. 2012) using 
the Harten & Yee TVD scheme as the convection difference operator in the fractional step 
method. The Wang et al. new high order finite difference method with subcell resolution 
procedure was developed for a single reaction to overcome the wrong propagation speed 
of discontinuities for stiff source terms. Here we apply the subcell resolution method on 
one of the reaction coefficients to observe its performance. 

Among the considered schemes, Fig. 2 indicates that the least dissipative scheme pre- 
dicts the shear and shock locations best when compared with the reference solution. The 
results indicate that TVDafi+split is slightly more accurate than WEN05-llf. This is due 
to the fact that TVDafi+split reduces the amount of numerical dissipation away from 
high gradient regions. The test of using the subcell resolution method by applying it to 
only one of the reactions in this multireaction flow, does not improve the performance 
over standard schemes. 

Figure 3 shows the shock velocity time dependence obtained on the four levels of grid 
refinement with Ax = 10 -3 to, 5 x 10 -4 m, 5 x 10 -5 to, and 2.5 x 10 -5 to. The shock 
velocity is computed by taking a numerical derivative of the shock location as a function 
of time with some smoothing. During the first 5 x 10 -6 sec the computed velocity has a 
strong dependence on the grid. After this initial duration, the computed velocities on the 
four grids asymptotically approach the same level. The result suggests that the major 
contribution to the error in shock location obtained on the coarse grid is due to the first 
5 x 1CV 6 sec. A similar behavior has been observed by other authors; see e.g., Jacobs 
(1994); Petrie-Repar (1997) for the perfect gas case. 

3.2. 2D EAST simulation results 

For the 2D case the computational domain in y is half of the 2D shock tube height. It 
has total length 8.5 to and height 0.0508 to (see Fig. 4). Other parameters and initial 
conditions of the high and low pressure regions are the same as for ID case. The bottom 
boundary is treated as an isothermal wall with the constant temperature T wa u = 300 K . 
The top boundary is treated as a symmetrical boundary condition. 
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Figure 5. 2D, 13 species EAST simulation using TVD for CFL = 0.7 on a grid with clustering 
between shear and shock. Temperature contours for tend, from top to bottom: 3.25 x 10 _5 sec, 
10~ 4 sec, 2 x 10 _4 sec. 


In the ^-direction both uniform and non-uniform grids are used for these simulations. 
For uniform grids the same strategy is applied in the ^-direction as for the ID case. For 
non-uniform grids the shock and shear locations are computed for each timestep, and the 
grid points are clustered in the cc-direction between the shear and shock locations with 
some tolerance to avoid interpolation errors during the re-grid process. This clustered 
portion of the grid moves as time advances. A grid stretching is also applied to smooth 
the transition from coarse to fine grid zone. At each time step the shock and shear 
locations are analyzed. If the shock/shear positions change by a prescribed distance, a 
re-grid is invoked and the data is interpolated onto the new grid. Note that in our case the 
interpolation is needed only in the vicinity of the transition between coarse and fine grid 
zones. The shock and shear locations are computed far away from the boundary layer. 
However, as the boundary layer develops, the shear layer becomes more curved so that 
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Figure 6. 2D, 13 species EAST simulation using TVD for CFL = 0.7 and t en d = 10~ 5 sec: 
Top Row: Three x-direction grid refinements: 601 x 121, 1201 x 121 and grid clustering between 
shear and shock in the x-direction (691 x 121 total). All y-grids use boundary grid stretching 
with a minimum of A y = 10 -5 m. Bottom Row: Two x-direction grid refinements: 1201 x 121 
and grid clustering between shear and shock in the x-direction (691 x 121 total). All y-grids use 
boundary grid stretching with a minimum of Ay = 5 x 10~ 6 m. 


the tolerance for the shear side of the grid-clustering zone is increased. In the y-direction 
all grids use the same stretching algorithm, so that the grid points are clustered in the 
vicinity of the boundary layer. 

The temperature contours computed for the times t en d = 3.25 x 10 -5 sec, t en d = 
10 -4 sec and t en d = 2 x 10 -4 sec are shown in Fig. 5. The solution is obtained on the 
grid with clustering between the shear and the shock. The grid spacing parameters are 
A X mirl = 5 X 10 -5 TO, A X m ax = 8 X 10 -4 TO, A y min = 10 -5 TO. 

For this 2D test case it is not practical to obtain a very accurate reference solution 
due to the CPU-intensive nature of the problem. Here, three levels of refinement are 
conducted. Figure 6 shows the computed temperature contour results at time t en d = 
10 -5 sec using CFL = 0.7 with TVD for three levels of x- and y-direction grid refinement. 
The top row shows three x-direction grid refinements of 601 x 121, 1201 x 121, and grid 
clustering between shear and shock in the x-direction of 691 x 121. The minimum grid 
step in the x-direction for the grid clustering is Ax m? ;„ = 5 x 10~ 5 to. All y grids use 
boundary grid stretching with a minimum of Ay = 10~ 5 to. The bottom row shows the 
same two x-direction grid refinements 1201 x 121 and grid clustering between shear and 
shock in the x-direction of 691 x 121 (fine block). All y grids use boundary grid stretching 
with a minimum of Ay = 5 x 10~ 6 to. Comparing the two rows of the grid refinement 
study indicates that refining the x-direction grid, keeping the y-direction the same, has a 



ID and 2D simulations for the NASA EAST Experiments 


9 



Figure 7. ID and 2D, 13 species EAST simulations using TVD for t en d = 3.25 x 10 5 sec: 
ID, Aa: = 5 x 10 -5 m (line 1), ID, Ax = 2.5 x 10~ 5 m, (line 2), 2D, 

A Xmin — 5 X 10 ” 771, Ay m i n - — 10 772 (line 3), 2D, A Xmin — 5 X 10 ” 772, Aymin — 5 X 10 772 
(line 4), 2D, A x m in = 2.5 x 10” 5 m,Ay m in = 10 -5 m (line 5). All ID simulations are on a 
uniform grid and all 2D simulations are on grids with clustering between shear and shock. 







T: 

2000 6000 10000 

14000 18000 



Figure 8. 2D, 13 species EAST simulation for t en d = 10 5 sec on the same grid with refinement 
between the shear and the shock: WEN05-llf for CFL = 0.2 (left), WEN05P-llf for CFL = 0.4 
(center) and TVD CFL = 0.7 (right). 


big effect on the locations of the shear and shock. This is due to the fact that aside from 
inside the boundary layer, the shear and shock are nearly one dimensional. However, 
comparing the last two columns of the grid refinement study indicates that refining 
the y-direction grid, keeping the ^-direction the same, has no effect on the locations 
of the shear/shock, but improves the boundary layer prediction. As in the ID EAST 
simulation, the discontinuity locations shift as the x-direction grid is refined, and the 
distance between the shear and the shock shrinks as the grid is refined. The shear and 
shock strengths are also different. Table 2 indicates the maximum shear and contact 
temperature for each set of grids. For the minimum grid stretching of Ay = 10~ 5 m, 
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Figure 9. 13 species ID EAST problem with zero source term. For notation see Figure 1. 


Grid N x 

601 

1201 

1201 

691 

691 

Cluster in x 

no 

no 

no 

yes 

yes 

Min Ay, m 

10 -5 

10“ 5 

5 x 10 -6 

10 -5 

5 x 10 

Shock 7 , K 

15,846 

18,851 

18,848 

25,098 

25,015 

Shear T rnax ,K 

11,301 

11,203 

11,203 

10,598 

10,598 


Table 2. Shock and Shear maximum temperature grid dependence at time t en d = 10 _6 sec. N x 
indicates the grid spacing in the ^-direction. The last two columns are grid clustering results for 
two different minimum j/-grid stretching. 


the maximum shear temperature is 11, 300 A', and the maximum shock temperature is 
15, 846AT for the 601 x 121 grid. However, the shear and shock strengths have a maximum 
shear temperature of 11, 200AT and a maximum shock temperature of 18, 851A' for the 
1201 x 121 grid. For the stretched grid the shear and shock strengths have a maximum 
shear temperature of 10, 600A' and maximum shock temperature of 25, 098A'. As we 
decrease the minimum grid stretching to Ay = 5 x 1CP 6 to, the shear and shock strengths 
have a maximum shear temperature of 11, 200 A" and a maximum shock temperature of 
18,848AT for the 1201 x 121 grid. For the stretched grid the shear and shock strengths 
have a maximum shear temperature of 10, 600 A' and a maximum shock temperature of 
25,015Ab Aside from the different shock/shear locations, the results in the last column 
show that the maximum temperature at the shock location is higher than the result 
indicated in the the middle and the first columns. 

To ensure that grid clustering with A x m i n = 5 x 1CP 5 to is enough in terms of shock 
location error, we performed a computation on a grid twice-refined in the ^-direction. 
The comparison of these results and a similar computation for the ID problem are shown 
in Fig. 7. The 2D result profiles are taken on the upper boundary of the computational 
domain, where they have the best correspondence with the ID case. The 2D results with 
the same Ax but different values of Ay almost coincide. The difference in the shock 
location for all five profiles does not exceed 1%. 

The comparison among three numerical schemes is shown in the Fig. 8. The three 
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methods are regular TVD for CFL = 0.7, regular WEN05 with Lax-Friedrichs flux 
(WEN05-llf) for CFL = 0.2 and a positivity preserving version of WEN05-llf (Zhang & 
Shu 2012), designated as WEN05P-llf, for CFL = 0.4. It is well known that WENO-llf 
is very diffusive. That results in larger errors in the shock location, as can be observed in 
Fig. 8. WEN05P-llf obtains results similar to regular WEN05-llf, but it appears to be 
stable even for higher CFL numbers, for which regular WEN05-llf obtains an oscillatory 
solution. 


4. Conclusions 

The present study demonstrates some important numerical challenges affecting the 
accuracy of ID and 2D numerical solutions in simulations of NASA EAST experiments. 
In the early stages of the time evolution, at around T= 10 -5 sec, the importance of 
obtaining very high accuracy of the discontinuities, in order to avoid an overestimation 
of the shock velocity, has been shown. For the 2D case, to obtain high resolution results, 
a moving grid with grid clustering only in key regions of the computational domain is 
needed in order to avoid the use of the full computational domain and unnecessary grid 
clustering for the entire time evolution. The cause of the observed grid-dependence of 
the numerical solution is not fully understood and requires further investigation. One 
conjecture for the spurious behavior might be the stiff source terms, or, at least, the grid 
dependence may be amplified when the stiffness of the considered Ll s is high enough. As 
discussed in (Wang et al. 2012; Yee et al. 2012) the level of grid- and scheme-dependence 
for obtaining the correct locations of discontinuities is normally dictated by the degree 
of stiffness of the source terms and the accuracy and amount of numerical dissipation 
contained in the scheme. Figure 9 shows the results of the same computation as in Fig. 1, 
but with zero source terms D s . A similar behavior is observed. In this case it appears that 
the problem is not physically realistic, especially when both the shear and shock jumps 
are extremely high. The present study also indicates the danger in practical numerical 
simulations for problems containing stiff source terms where there is no reliable means 
of assessing the accuracy of the computed result other than by extreme grid refinement, 
which may be beyond the capability of current supercomputers. Another alternative 
would be to develop methods that obtain the correct speed of discontinuities on coarse 
grids, e.g. using ideas similar to (Wang et al. 2012). This approach might be very useful 
when there is a need to perform 3D computations for such experimental facilities as 
EAST. Future investigation is planned. 
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